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We investigate the transport properties of a one-dimensional superconductor-normal metal- 
superconductor (S-N-S) system described within the tight-binding approximation. We compute 
the equilibrium dc Josephson current and the time-dependent oscillating current generated after the 
switch-on of a constant bias. In the first case an exact embedding procedure to calculate the Nambu- 
Gorkov Keldysh Green's function is employed and used to derive the continuum and bound states 
contributions to the dc current. A general formalism to obtain the Andreev bound states (ABS) of 
a normal chain connected to superconducting leads is also presented. We identify a regime in which 
all Josephson current is carried by the ABS and obtain an analytic formula for the current-phase 
relation in the limit of long chains. In the latter case the condition for perfect Andreev reflections 
is expressed in terms of the microscopic parameters of the model, showing a limitation of the so 
called wide-band-limit (WBL) approximation. When a finite bias is applied to the S-N-S junc- 
tion we compute the exact time-evolution of the system by solving numerically the time-dependent 
Bogoliubov-deGennes equations. We provide a microscopic description of the electron dynamics not 
only inside the normal region but also in the superconductors, thus gaining more information with 
respect to WBL-based approaches. Our scheme allows us to study the ac regime as well as the 
transient dynamics whose characteristic time-scale is dictated by the velocity of multiple Andreev 
reflections. 

PACS numbers: 



I. INTRODUCTION 

In the last few years nanoscopic Josephson junctions 
have been widely studied both theoretically and exper- 
imentally as possible candidates to provi de an alt erna- 
tive technology to silicon-based electronic d 1 * 2 * 3 * 4 * 5 -^. Spe- 
cial attention has been paid to the analysis of super- 
conducting atomic-size quantum point contacts (SQPC^l 
like single-level quantum dots and nanowires. Among 
the most striking features experimentally observed we 
mention the subgap structure in the current-voltage 
characteristics driven by multiple Andreev reflections^, 
the single-electron tunneling through discrete electronic 
states^, and the nanoscopic dc Josephson current!^. 

Within the so-called Hamiltonian approac}^^ it is pos- 
sible to provide an accurate microscopic description of 
these systems, where some relevant length scales (Fermi 
length, size of the junction, etc.) are comparable. This 
approach relics on tight-binding-likc Hamiltonians and 
has the advantage to treat the tunneling Hamiltonian 
describing the SQPC to all ordersPHH. In SQPC the 
Andreev bound states (ABS ^21111 pl a y an important role 
since they can carr y an important amount of dc Joseph- 
son currentF- 4 * 15 * 16 [ Such states origin from multiple An- 
dreev reflections occurring at the superconductor-device 
contact and come in pairs, one above and one below the 
Fermi level, carrying opposite supercurrents. In spite of 
the large theoretical effort in studying the dc Joseph- 
son regime in SQPC, a proper description of extended 
junctions is still lacking since the electrodes degrees of 



freedom have been so far absorbed in an approximate 
frequency-independent pairing and on-site potentials at 
the boundaries of the central regiorP^ 18 l 19 [ 

The calculation of the ac Josephson current is more 
involved. At present the ac regime has been studied us- 
ing Floquet-based methods co mbined w ith nonequilib- 
rium Green's function tectinique d 11 * 20 * 21 ^. This approach, 
however, is limited to the dc bias case and other interest- 
ing time-dependent driving fields, like ac bias or voltage 
pulses, cannot be addressed. A possible alternative ap- 
proach is the one based on the real time-propagation but, 
so far, only normal metal-quantum dot-superconductor 
junctions have been studied 22 . 

In this paper we investigate the transport properties 
of a one-dimensional (ID) superconductor-normal metal- 
superconductor (S-N-S) systenpH composed by a normal 
tight-binding chain embedded between two ID supercon- 
ductors described by the Bogoliubov-deGennes Hamilto- 
nian. We will study both the static dc Josephson current 
J and the time-dependent oscillating current generated 
after the switch-on of a constant bias. In the dc case we 
employ an exact embedding procedure and calculate the 
three different contributions to J, carried by the ABS, 
the normal bound states (with energy below the bottom 
of the band), and the continuum states. We show that if 
the pairing potential is larger than half the bandwidth of 
the normal region, all Josephson current is carried by the 
ABS's. In this regime we are able to extend the results 
by Affleck et alP^in the limit of long normal region. The 
use of the exact embedding self-energy allows us to re- 
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late the phenomenological paring potential of ReflTTlwitli 
the microscopic parameters of the model, thus obtaining 
a condition for perfect Andreev reflections in term of the 
physical order parameter A. In addition we highlight a 
limitation of the commonly used wide-band-limit (WBL) 
approximation . 

When a finite bias is applied to the S-N-S junc- 
tion, we compute the exact time-evolution of the system 
by solving numeric ally the time-dependent Bogoliubov- 
deGennes equation d 13 * 24 * 25 l This is done within the so- 
called partition-free approach, in which the S-N-S system 
is assumed to be contacted a nd in equilibrium before the 
external bias is switched onJ 26 * 27 *. Explicit calculations 
are performed in the case of superconducting leads of fi- 
nite length. However, as already discussed in Refl2"5l the 
electrodes are long enough to reproduce the time evo- 
lution of the infinite-leads system. The above approach 
gives us the possibility to explore the transient dynamics 
and provides a time-dependent picture of the Andreev re- 
flections. In the long-time limit we recover the expected 
oscillating current, whose Fourier transform displays con- 
tributions from different harmonics of the fundamental 
Josephson frequency. By extracting the dc component of 
the oscillating current, we are also able to reproduce the 
subgap structure in the current- voltage characteristics. 

The paper is organized as follows. In Section |TT] we 
introduce the model Hamiltonian and briefly recall the 
Nambu and Bogoliubov-deGennes formalisms. In Sec- 
tion |III| the equilibrium Josephson current is studied by 
means of an exact embedding procedure. Numerical re- 
sults for short junctions are reported in Section [TV| while 
the limit of long normal regions is analytically carried 
out in Section [V] In Section [Vl| we investigate the time- 
dependent regime. Two Appendices corroborate the ana- 
lytic derivations. Finally summary and main conclusions 
are drawn in Section IVIII 



The Hamiltonians for the Left/Right (L/R) supercon- 
ducting leads has the general form 
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6 XQ ^}aJ. C <?aT 



, (3) 



where A Q is the pairing potential in lead a — L, R with 
corresponding phase Xa- The one-particle energies e q 
span the range (— W, W) where 2W is the lead band- 
width. We assume that the tunneling between the su- 
perconductors and the normal region occurs only via the 
boundary sites of the chain and model Ht as 



H T = ^V q \c\ L]Cl] + c\ R] c M ] 



h.c. 



(4) 



In the last term of Eq. ([I]) /i is the chemical potential and 
N-t/i is the number of electron/holes with spin f / J,. 
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FIG. 1: Scheme of the S-N-S junction. For illustration the su- 
perconducting leads are ID chain with nearest-neighbor hop- 
ping ts [i.e. e q — 2tscosg in Eq.Q] and on-site pairing 
potentials Al = Ar — A with \l = Xfl — 0. The hopping 
integral between the boundary sites of the superconducting 
and normal regions is tr [i.e. V q = tT sin in Eq. (J4J , 

where A is the number of sites in the leads] . 



II. THE MODEL 

We consider a hybrid S-N-S system consisting of a nor- 
mal region contacted to two superconductors, as illus- 
trated in Figjl] In the Bogoliubov-deGennes formalism 
the annihilation (creation) fermion operators anni- 
hilates (creates) electrons of spin up, while the annihi- 
lation (creation) fermion operators cj^ annihilates (cre- 
ates) holes of spin down. In order to avoid confusion we 
put a tilde on the hole-operators. The Hamiltonian of 
the system is described by 

H = H N + H L + H R + H T -n{N 1 -N l ). (1) 

In this work we consider normal regions consisting of a 
tight-binding chain of length M and nearest neighbor 
hopping tjv with Hamiltonian 



M-l 



H N — t 



N 



E K 



(2) 



The time-dependent current! 2 ^ at the a = L, R inter- 
face is 

I a (t) = 2j2V q ReTr[G< qa (t,t)] , (5) 



where the Nambu lesser Green's function is defined as 

[G < (t 1 ,t 2 )] m ,n = G^ ln (ti,t2) 



lT (ti)c„ T (i 2 )) {$ ml {ti)c n] {t 2 )) 



A {ti)c ni {t 2 )) {^ mi {ti)c ni {t 2 )) 



(G) 



In the above definition the indices to, n denote either a 
site in the normal chain or a g-state in the a = L, R 
lead. We observe that the off-diagonal components of 
the Green's function can be interpreted as spin-flip prop- 
agators in the effective Bogoliubov-deGennes space. The 
retarded, advanced and greater Green's functions are de- 
fined in a similar way as in Eq. ^ . 

The rest of the paper is devoted to the calculation of 
I a (t). First we will focus on the equilibrium problem 



3 



and calculate the dc Josephson current J — -Zl(O) = 
Ir(0). Then we apply a finite bias voltage across the 
junction and compute numerically the time-dependent 
current /l(£) at the left interface. 



III. DC JOSEPHSON CURRENT 

The dc Josephson current J = iz,(0) is obtained from 
Eq.Q with an equilibrium lesser Green's function and 
reads 



J = 2Re 



dio r 
2tt" 



-Tr 



L q 



(7) 



In equilibrium the lesser Green's function is related to the 
retarded/advanced Green's function via the fluctuation- 
dissipation theorem 



G<(u,) = -/H[G»-G» 



(8) 



where / is the Fermi distribution function. In the follow- 
ing we work at zero temperature. This means that the 
effective pairing potential in Eq.Q corresponds to the 
BCS gap at T = 0. The entire formalism remains valid 
at finite temperature T, provided that the order parame- 
ter A corresponds to the BCS gap at T. The dependence 
on temperature of the current J is mainly due to the 
change of A(T), since the Fermi function / remains close 
to a theta function for T < A. This is supported by 
the results shown in Fig(6]which agree well with previous 
studies on temperature dependence of J"™. 

By exploiting the Dyson equation for the re- 
tarded/advanced Green's function the Josephson current 
J can be expressed in terms of the embedding self-energy 



j a 



r/a 



(9) 



as 



J = 2Re / ^tP* { [Sl,i( u ) - Gj,i(w) a z } 

. (io) 

where cr, is the third Pauli matrix and q T ' a is the Green's 

z ^-aq 



function of the isolated a lead. We observe that Eq.(10) 



is valid for any S-N-S system provided that the S-N hop- 
ping occurs only at the two boundary sites of the normal 
region. The general expression for the embedding self 
energy is 



S 7 » = 



m a (uj±iri) A„(w ± ir\)e %Xo 
A a (uj ± ir])e~~ lXa m a {uj±irj) 



(11) 

where m and A are the effective on-site and pairing po- 
tentials. In the case of ID superconducting leads with A 
sites (see Fig{T| 



2ts cos q 



V q = t T \j -sm ,j 



(12) 



and the self-energy at — is (see Appendix [A} 



.ID 



t\ y/A 2 a - Z 2 _ - Z2 



4t| 



'2i| 



Ai D (z) = A a 



4*| 



2t% 



(13) 
(14) 



where z is a complex frequency and the infinite A limit 
has been taken. The WBL result is easily recovered 
by defining the tunneling rate T = 2t^/ts, expanding 
Eqs.( fl3|5 l in powers of z/ts and A/ts and retaining 
only the zero-th order term. In this way one gets 



A WBL (z) 



r z 
~2jEi^ 
r a q 

2 7AT^ 



which for z = uj + iij yields the commonly used WBL self- 
energy 2 WBL . We would like to observe that evaluating 
E WBL at the Fermi energy u> = /j, = one finds that 
m WBL _ q and tllat the p a i r i n g potential A WBL oc T 

is independent of the order parameter A . In Section [V] 
we will discuss the implications of this feature for long 
normal chains. 

In the rest of the Section we do not assume any spe- 
cific form of the embedding self-energy and present a 
general procedure to calculate the dc Josephson cur- 
rent of Eq. ( 10 1 . For practical purposes we split the in- 
tegral in Eq.(10) in three different energy regions and 



identify the contributions of the normal bound states, 
Andreev bound states and continuum states (see Fig|2|. 
The energy range is (-^W 2 + A^ ax , -A min ) for the 




~Y Y 

normal bound states band continuum Andreev bound states 

FIG. 2: Schematic representation of the density of states of 
the S-N-S system. The three spectral regions corresponding 
to normal bound states, Andreev bound states and continuum 
states are displayed assuming Al ~ Ar — A. 



filled continuum states, 



A 2 

max 



(-A min ,0) for the filled ABS's 
) for the filled normal bound 
max{A L ,A fi } and A min = 



and (—oo, — \/W 2 
states, where A n 
min{Ai, A^}. Thus letting j(u>) be the integrand func 
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tion in Eq.(lO) the total Josephson current reads 



J 

J cont 
J abs 
J nbs 



= J< 



cont "1" J, abs ~h J nbs 



-A„ 



doj 



2tt 



-y/W 2 +Al 



2^ 



(15) 
(16) 

(17) 

(18) 



The nature of the above decomposition is illustrated in 
Figj3] where the the integrand function j{uj) is displayed 
for a ID S-N-S junction at a fixed value of \ — Xl — 
Xr = 7r/3. In Figj3]we have chosen the superconducting 
gap A about one order of magnitude smaller than the 
leads bandwidth W — 4tg. This is done in order to 
highlight the contribution coming from the normal bound 
states, although we expect that it becomes less and less 
important as A/W — ► 0. We would like to emphasize, 
however, that the normal bound states play a crucial role 
in a self-consistent calculation of the total current, since 
the effective mean-field potentials depend on the density 
and in the central region the contribution of the normal 
bound states is certainly not negligible. Indeed the total 
number of particles iVjv in the central normal region is 



N N = -i 



r j M 



(19) 



where the integrand function in the above equation has 
a similar structure as in Fig|3j 




-0.5 



FIG. 3: Integrand function j(uj) in Eq.p5|| for a ID S-N-S 
junction with M = 4, ts = tr = 1, ijv = 1.2, Al = Ar = 
0.6, x — Xl — Xr = t/3. The two dashed vertical lines 
correspond to w = —\/W 2 + A^ ax and ui — — A m j n and mark 
the boundaries of the three integration regions. A broadening 
77 = 10 -4 has been used to give to j a finite width around the 
bound states. Energies are in units of \ts\- 

It is worth noticing that the function j(u>) is propor- 
tional to a Dirac delta around the bound states. There- 



fore the numerical integrals in Eqs.(17l and (18 1 must 
be computed with care. An efficient alternative way to 



calculate with high accuracy J a b s and J n b s consists in 
realizing that 



J, 



abs 



gV ^abs 

4- dX 



J, 



nbs 



^ d X 



(20) 



where E^ s and are the energies of the filled An 



dreev and normal bound states respectively. Eq.(20l 



follows directly from the Hellmann-Fcynmann theorem, 
which in this case can be exploited since E^ s and E^ s 
are the eigenenergies of the Hamiltonian in Eq. ([lj . By 
a simple gauge transformation the phase \ can be trans- 
ferred to the hopping integrals V q in Eq.Q, and hence 



the derivative of E^ s and E^' s with respect to x yields 
the average of the current operator over the Andreev and 
bound eigenstates. In the following we derive and elegant 
formula to calculate E^ and E^. 

The bound state energies can be obtained by solving 
the self-consistent 2M x 2M secular problem 

reffCE-M./. \ _ TPtJ. \ (21) 



(n) 



H°«{E)\iP E ) = E\^ E ) 



with \E\ < A min (Andreev) and \E\ > \/W 2 + 
(normal), and 



Hf(E) = H 



N 



+ m R (E)[ 

c mt Cm T + c\ n CMi] 



Ml M \ 1 ^Ml K 

MI^mi t s IXr c'mi c mi\ 



A R (E) [e^c\ ~ CMi + e -**A 



+ m L (E) [c f 1T Ci T +4;Ci4.] 

+ A L (E) [e^c^cu + e-^c^c 1T ], (22) 



with m and A as in Eq.(ll). In the effective Hamil 



tonian the on-site and pairing potentials at the bound- 
ary sites 1 and M are renormalized by the embedding 
procedure. In order to simplify the algebra we define 
the momenta k such that E = 2tx cos(k) and assume 
Al = Aji = A max = A m ; n = A, which also implies 
A L (E) = A R {E) = A fe and m L (E) = m R (E) = m k . 
In Appendix [5] we describe in detail how the eigenvalue 
problem in Eq.(|21[) is analytically solved to yield the fol- 
lowing equation for the momenta k 



= 4sin 2 fc(M+ 1) + ( 



-l) M 2t 2 N Al cos x sin 2 k 



~ 2t N A 2 k sm 2 (kM) 
~ m 2 k [2t N sm(kM) - 



- Al sin 2 k(M 
nik sin k(M — 



-1) 
I)] 2 . 



(23) 



where x — Xl — Xr- The bound state energies are found 
by solving Eq.(23l and retaining only the values of k 
for which |2tjvCOsfc| < A and \2t N cosk\ > \/W 2 + A 2 . 
We observe that in general the variable k is complex, see 
Appendix |B| 

We would like to end this Section by commenting two 
limiting cases. For an isolated normal region (m k — 
Afe = 0), the allowed momenta are simply k = irj/(M + 
1), j = 1,...M, as expected. We also observe that if 
we set TTife = and assume a constant pairing poten- 
tial Afe = A, the above equation reduces to Eq.(3.3) of 
RefJTTl 
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IV. NUMERICAL RESULTS FOR J( x ) 

By following the approach described in the previous 
Section, we specialize to the case of half-filled ID leads 
as in FigjT] (to = to 1d and A = A 1D ) and numerically 
evaluate the dc Josephson current. In Fig|4]we show the 
current J as a function of \ as well as the three dif- 
ferent contributions J cont , J abs , J nbs for a chain of M = 8 
sites. We notice that there is an optimum value of fjv [see 




FIG. 4: Total Josephson current J (solid curve), J CO nt (dotted- 
dashed curve), J a bs (dashed curve) and J n bs (dotted curve) as 
a function of x = XL — XR f° r a S-N-S junction with M = 
8, ts = tr = 1, Ax, = Ah = 0.6. The panels (a) to (f) 
correspond to tjv — 1.5,1.2,1.0,0.744,0.6,0.3. Energies are 
in units of \ts\- 

Figj4] panel (d)] at which there is a non trivial cancella- 
tion of the non-linear contributions J CO nt and J a b s and 
J(x) becomes a straight line. In this regime the Joseph- 
son current is also maximized for every value of \. We 
have further investigated this instance and found that for 
any given A = A^ = Ar, there exists an optimum value 
of tj\[ — tjf at which this property is observed. In the 
left panel of Fig.([5| we plot as a function of A for 
the same parameters as in Fig(4] We have also observed 
that tjf is quite insensitive to the size M of the normal 
region. In the right panel of Fig. (JEJ we display the cor- 
responding critical current j( b \ir), i.e. the value of the 
Josephson current reached at x — 7r - The linearization 
of the current-phase relation is also known as the Ishii's 
sawtooth behavior^ and corresponds to perfect Andreev 
reflection^. We notice that the ABS contribution sat- 
urates the dc Josephson current for % = 0.3, see Fig|4] 
panel (f). In the next Section we consider the limit of 



0.14 
0.12 
_ 0.1 

f o.os 

0.06 
0.04 



FIG. 5: 



(Left panel) tf ] 



as a function of A = Al = Aj 



(Right panel) Critical current J(tt) as a function of A cal- 



culated at tjv 



The rest of parameters are M 



ts = tr = 1. Energies are in units of \ts\- 



long chains and identify a regime for the occurrence of 
this saturation. 

Finally we show that within our approach it is possi- 
ble to reproduce the crossover of the current-phase rela- 
tion between short and long S-N-S junctions, already dis- 
cussed in previous workiPS. In Fig|(3]we display J(x) both 
for a short junction (M — 1) as well as for a long junc- 
tion (M = 51). It appears that for large superconducting 
order parameter the current-phase relation evolves from 
a sin(x/2)-shaped curve^ to a straight line by passing 
from to M = 1 to M = 51. These results are in agree- 
ment with the findings of Ref.30 where the change of the 
order parameter A is due to a change of temperature. 



JOc) 



(a) 
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7T 




n/2 x " 

FIG. 6: Total Josephson current J as a function of \ = XL — 
X r for a short junction with M = 1 (panel a) and for a 
long junction with M = 51 (panel b) for different values of 
A i, = A R = A = 0.1,0.01,0.005,0.001 (from top to bottom). 
The rest of parameters are ts = 1, tjv = 3.8, tr = 2. For 
clarity, the curves corresponding to A = 0.01 and A = 0.005 
have been multiplied by a factor 5, while curves corresponding 
to A = 0.001 by a factor 10. Energies are in units of |is|- 



V. LIMIT OF LONG NORMAL REGION 

In this Section we study the Josephson current in the 
limit of long chains. By numerical inspection we have 
verified that J = J a b s for t/v < A/2, i.e. all the Joseph- 
son current is carried by the ABS's. In this regime the 
number of occupied ABS's equals exactly the number of 
sites M of the tight-binding chain. Thus the ABS's con- 
stitute a local basis set with a good approximation. As 
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a consequence no normal bound states occur, while the 
amplitude of the current carrying continuum states is ex- 
ponentially suppressed in the normal region. The current 
J = J a bs is obtained by calculating the contribution E^ s 
of the ABS's to the total energy 



Ax) 



dx 



(24) 



(k) 

To calculate the energy E ah ' s = 2t/v cos A: of a single ABS 
it is convenient to write 



7TJ 



M + 1 M + 1 



J 



1,...M. 



(25) 



where 5k is a fc-dependent phase-shift. Following Refs,17, 
[35] the total ABS energy can be expressed as 

71 dk 



E™ = 2(21/ + 1) / 

JO 





dh a.bs 


/O 7T 





1 7TWi? 

2M + 1 



[4,+ + S k ,-]f(E^ s ) 



Sk F ,+ \ , ( far,- 



(26) 



with A fcF = A^. We would like to stress that Eq.(30l 
has been obtained starting from a microscopic model 
Hamiltonian, i.e. without resorting to phcnomcnologi- 
cal effective on-site and pairing potentials. The relation 
between the effective pairing potential A kF and the mi- 
croscopic order parameter A in Eq. ( 14 1 allows us to 



discuss some relevant limiting cases in terms of physical 
quantities. 

In Fig. ^ we plot the Josephson current in Eq.(29l 
using for the phase-shift the result in Eq.(30). We fix 



the values of the hopping parameters to be t^ — 0.618, 
ts = It — 1 and study how the current-phase relation 
depends on A. We notice that for A = 1 the current is 
linear in the ranges [0, 7r) and (tt,2tt], with a sharp dis- 
continuity at x = 7T- This is the Ishii sawtooth behavior^ 
already mentioned in the previous Section. In that case, 
however, the sawtooth behavior was the result of a per- 
fect cancellation between the contribution of the contin- 
uous states and of the ABS's. We also verified that the 
Josephson current calculated by means of the brute-force 
numerical evaluation of Eq.(10) at < A/2 is in excel- 
lent agreement with the current evaluated as in Eq. ( 29 1 
already for M > 10. 



where 5 k ,± correspond to the two branches of ABS's, 
vf = 2t n sin kp is the Fermi velocity and kF is the Fermi 
momentum. For large M the momentum k is a contin- 
uous variable in the range (0, tt) and the phase-shifts 5k 
can be determined by inserting Eq. ( 25 1 in Eq. ( 23 1 and 
expanding in powers of X/M. To lowest order Eq.(23) 
reduces to 



(a k sin 5 k 



where a k 
b k ■■ 
2(t N 



f2 

>'N 



b k cos 5 k ) 
- A 2 )cos(2fc) 



c k . 



(27) 

Itjssrrik cos k, 



= 2tjqmkSm.k — {jn\ 

AfcSin/« 2 (l — cosx). The solutions of Eq.(27l read 



A 2 .) sin(2fc) and c 



c . b k 1 

dk ± = — arctan — ± - arccos 
a k 2 



1 - 



(28) 



Eq.(28 I provides a generalization to nonvanishing on-site 
potential of the phase-shifts found by Affleck et alP^. 
Inserting Eq.(28) in Eq.(26l the Josephson current is ob- 
tained from Eq.p4]). We notice that the combination 
5k. + + 5k,- is independent of the phase difference x f°r 



any A fc and m k . 
reads 



Therefore the dc Josephson current 



Ax) 



ttvf d 
M+ldx 



>k F ,- 



(29) 



Below we specialize the analysis to ID leads at half-filling 
(kp = tt/2). In this case m k u = 0, see Eq.(13|, and one 
can show that 



fik f 



±i arccos 

-j 



(cosx - 1) 



(t 2 



N 



A? 



(30) 




n/2 
X 



FIG. 7: Josephson current as in Eq.(29l for different valu es 
of A for t N = 0.618, ts = t T = 1. For A = 1 Eq.pi) is 
fulfilled. J(x) is in unit of vf/{M + 1). Energies are in units 
of |t s | . 



As shown in ReffTTl the linear behavior of J is due to 
perfect Andreev reflections which occur for ijy 
i.e., 



A 2 -A). 



(31) 



We recall that the above current corresponds to the total 
Josephson current only for tN < A/2, which, together 
with Eq.(31l, implies 



t N < 



(32) 



Equations ( 31|32 1 establish a regime in which the Joseph- 
son current is entirely carried by the Andreev bound 
states via perfect Andreev reflections. 
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Before concluding this Section we would like to observe 
that in the WBL approximation the condition for perfect 
Andreev reflection implies t^ = A^ BL = T/2, which 
does not depend on the order parameter A. The same 
limitation of the WBL approximation emerges in the cal- 
culation of the phase-shifts, see Eq.(30). Therefore the 



use of WBL self-energies in superconducting transport 
through long normal chains does not allow to study the 
dependence of the current-phase relation on the physical 
order parameter. 



VI. AC JOSEPHSON CURRENT 

In this Section we consider the time-dependent current 
flowing through the S-N-S junction after the switch-on 
of a dc bias voltage. In order to get a sensible transient 
regime, we adopt the so-called partition-free approach, 
in which the S-N-S system is assumed to be contacted 
a nd in equilibrium before the external bias is switched 
OI j26 | 27 | rp ne numer i ca i results contained in this Section 
are obtained by computing the exact time-evolution of 
the system described in Eq.Q with finite ID supercon- 
ducting leads of length A (see Fig[T]). Without loss of 
generality we switch on the bias at t = 0. The biased 
Hamiltonian at positive times reads 

H{t) = H N + H L (<) + H R {t) + H T - M ( A T - JVj.) , (33) 

where 



+ A p i(x*+2U a t) t ? 



A e -i(X a +2U a t)~i 



qOt] 



(34) 



and U a are the dc bias voltages applied to lead a. We 
denote with H(0) the matrix representing the equilib- 
rium Hamiltonian H of Eq. Q projected over one-particle 
states and with H(t ) the corresponding matrix represent- 
ing H{t) of Eq.(34) for t > 0. The generic element of 



H(t > 0) is a 2 x 2 matrix in the Bogoliubov-deGennes 
space 



[H(i)] m: „ — 



H„ 



-H„ 



.(«) 

,n(t) 



(35) 



where m,n — 1, 2 A + M. According to the partition- 
free approach, we first calculate the equilibrium config- 
uration of the contacted system by solving the secular 
problem 



£[H(0)] r 



u k (n) 
v k (n) 



E {k) 



u k {m) 
v k (m) 



(36) 



and construct the initial lesser Green's function 



[G < (0,0)] m , n = 



i[/(H(0))] m ,„ 

u* k [m)u k (n) u* k (m)v k (n) 
v* k {m)u k (n) vl{m)v k (n) 



.(37) 



The initial states are then propagated in time according 
to the time-dependent Bogoliubov-deGennes equations 



i—u k (m,t) 



[H m . n (t)u k (n, t) + A m>n (t)v k (n, t)} 

n 

i^rv k {m,t) = y~l \-H m<n {t)v k {n, t) + A min (t)u k (n, t)] , 

n 

(38) 

which are solved by 

: ^[ re -ij?*-BW 



u k (m,t) 
v k (m,t) 



u k (n,0) 
v k (n,0) 



(39) 

with initial condition Ufc(m,0) = u k (m) and v k (m, 0) = 
v k (m) and T the time-ordering operator. The lesser 
Green's function G < (t 7 t) has the same form as the r.h.s 
of Eq.(37l with u k (m) and v k (m) replaced by u k (m,t) 
and v k (m, t) . Expressing the time-dependent wavefunc- 
tions as in Eq.d39| it is straightforward to show that 



G<(t,t) =Te- % ti dT ^ {T) G<(0,0) Te l ^ dTS 



(t) 



(40) 



We notice from Eq.(34) that H(t) has an explicit time- 



dependence (the time-dependent phase of the order pa- 
rameter) and hence the evolution operator is not the ex- 
ponential of a matrix albeit the bias is constant in time. 
This problem is solved by discretizing the time and calcu- 
lating the evolution of the lesser Green's function within 
a time-stepping procedure 

G^tj,^) « e' im ^ )st G^-i.tj'-i) e^ )5t , (41) 

where tj = jdt, St is a small time step and j a positive 
integer. The time dependent current at the left interface 
is calculated from Eq.(|5]). The above approach allows us 
to reproduce the time evolution of the infinite-leads sys- 
tem up to a time T max ps 2A/v, where v is the maximum 
velocity for an occupied one-particle state. For t > T max 
high- velocity particles have time to propagate till the far 
boundary of the leads and back, yielding undesired finite- 
size effects in the calculated currenPSJ. For this reason 
we set A such that 2A/v is much larger than the time at 
which the stationary oscillatory state is reached. 

In FigjS] we plot the time-dependent current through 
a single-dot junction (M — 1) for different values of the 
superconducting order parameter Al = Ar = A, rang- 
ing from to 1. In panel (b) we display a magnification 
of the transient regime. It appears that the transient dy- 
namics becomes slower as A is increased. This is due to 
the fact that at bias U ~ 2A/n, an incident electron com- 
ing from the left superconducting lead undergoes about 
n Andreev reflections inside the central region before be- 
ing transmitted to the right lead. We also verified that, 
at fixed A, the transient timescale grows by reducing the 
bias voltage (not shown). A qualitatively similar behav- 
ior is observed in Fig(9j where the hopping in the super- 
conducting leads is taken about two orders of magnitude 
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FIG. 8: Current I hit) through the left interface for different 
values of = A# = A = (thin solid curve), 0.1 (dotted 
curve), 0.5 (dotted-dashed curve), 0.7 (dashed curve), 1 (thick 
solid curve). The rest of parameters are M = 1, A = 80, 
St = 0.2 t s = tr = -1, XL = XR = 0, U L = -U R = 0.25. 
Panel (b) displays a magnification of the transient regime for 
< t < 15. Energies are in units of \ts\, while time and St 
are in units of l/|ts|. 



larger that all the other energy scales, in the spirit of the 
WBL approximation. Another interesting observed fea- 
ture is that the dc component of the current II in Fig{8| 
displays a non-linear behavior with A. In particular II 
increases with A passing from to 0.5, but decreases by 
further increasing A form 0.5 to 1. Such behavior, how- 
ever, is not seen in Fig|9] where II is a monotonically 
decreasing function of A. 

At long time the current 1^ it) displays the well known 
ac Josephson behavior, with persistent oscillations at 
multiple frequencies of the fundamental Josephson fre- 
quency Uj = 2(Ul — Ur). To investigate the stationary 
oscillations we performed a discrete Fourier transform 
of in the time window (T m i n , T max ) where T mm is 

much larger than the transient timescalc. Denoting with 
Nf the number of time steps in the time window . the 
Fourier components of are defined according td 28 * 33 ! 



Nf 



N f — 

J ri — — 



(42) 



where u„ = 2irn/(Nf6t). In Fig 10 we plot the dissipa- 
tive contribution In{to n ) = 2Rc/(oj„) and the nondissipa- 



FIG. 9: Current Ih{t) through the left interface for different 
values of A^ = A^ = A = (thick solid curve), 0.25 (thin 
solid curve), 0.35 (dotted curve), 0.5 (dotted-dashed curve), 
0.75 (dashed solid curve). The rest of parameters are M = 1, 
A = 6000, St = 0.1 t s = 100, t T = 4.47 (i.e. T = 2t%/t s = 
0-4), XL = XR = 0, U L = -Ur = 0.5. Panel (b) displays 
a magnification of the transient regime for < t < 6. The 
time propagation has been obtained by retaining only the one- 
particle states in Eq.( |34[ ) with energy —10 < e q < 10. We 
have checked that within this choice the results with A = 
perfectly agree with ones of Ref 1271 obtained within the WBL 
approximation. Energies are in units such that T = 0.4, while 
time and St are in units of 1/r. 



tive one Jnd(wh) = — 2Im/(w„) to the currentP^S. The 
first four harmonics are clearly visible and the fundamen- 
tal component is the dominant one. We also observe that 
the amplitude of the harmonics is not a monotonically de- 
creasing function of the frequency. The above procedure 
provides an alternative method to perform the spectral 
decomposition of the ac Josephson current. Our time- 
dependent approach is not limited to dc biases and the 
same computational effort is required to study ac or more 
complicated time-dependent biases. From our numerical 
time-dependent simulations, it is also possible to extract 
the cur rent- voltage characteristics of the junction. In 
Fig.(ll I we show 1^ as a function of the applied dc bias 
for a S-S junction (M = 0). The system consists of two 
ID superconductors connected to each other via a hop- 
ping integral tr between the boundary sites of the L and 
R leads. We observe a well defined sub-gap structure 
characterized by current kinks at Ul — Ur — 2A/n, a 
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0.02 



-0.02 
-0.04 
-0.06 
-0.08 



(a) 




FIG. 10: Non-dissipative coefficient Ind (panel a) and dissipa- 
tive coefficient Io (panel b) obtained from the discrete Fourier 
transform of Ih{p) as described in the main text. They are 
calculated using 2000 equidistant points of II it) — II with t 
in the range (50, 140). In this plot A = 150, St = 0.05, A = 1 
and the Josephson frequency is ujj = 2(Ul — Ur) = 1. The 
rest of the parameters are the same as in Fig[8] Energies a 
and frequency are in units of \ts\- 



0.15 




0.05 



FIG. 11: Current-voltage characteristics (Tl vs Ul) of the S-S 
junction for different values of the hopping tr- The rest of 
parameters are Ur = 0, ts = —1 and A = 0.5. The vertical 
dotted lines denote the values Ul = 2A/n (n=l,2,3,4) at 
which multiple Andreev reflections are expected. Energies 
are in units of \ts\- 



feature already poin ted out in previous works within the 
WBL approximatiorfSEM]] \y e nave a \ so checked that 

if the WBL is modelled with ID leads (i.e. by taking 
t s > 1 and t T = y/Tts/2 with finite L), we numeri- 
cally recover the current- voltag e characteristics already 
obtained in previous work£LU22l. 

Finally we have computed the time-dependent evolu- 
tion of the spin-up electron density according to 



(t) = -i([G<(t,t)] r 



(43) 



where m denotes a site of the S-N-S system and the ma- 
trix element (. . . )i,i is taken over the Nambu space. We 
stress that our approach allows us to determine n m ^(t) 
not only in the normal region, but also inside the su- 
perconducting leadiP^. This is a clear advantage with 
respect to the WBL approximation, in which only the dy- 
namics of the normal region can be described. In Fig |12| 
we show the density variation <5n m f it) — n(t) m -\ — n m j (0) 
as a function of the atomic position m along the ID S- 
N-S system and time. In this case a long junction with 




40R 



FIG. 12: Contour plot of the time-dependent variation of den- 
sity for spin-up electrons 5n m |(t) = n(t) m i ~ n(0)mf as a 
function of the atomic position m along the ID S-N-S system 
(x axis) and time (y axis). Srn-\(t) is displayed for the first 
40 sites in both leads and inside the M = 20 sites of the nor- 
mal region. The rest of parameters are A = 100, St — 0.3, 
t s = 1.2, t N = l,t T = 0.8, A L = A R = 0.2, X L = XR = 0, 
Ul = 0.3, Ur = 0. Energies are in units of |tjv|, while time 
and St are in units of l/|ijv|. 




20R 40R 



FIG. 13: Same as Figjl2] The model parameters are: M = 21, 
A = 200, St = 0.3, t s = l,t N = 1, tr = 1-104, A L = A R = 
0.4, xl = XR — 0, Ul = 0.2, Ur = 0. Energies are in units of 
| £jv 1 , while time and St are in units of l/|ijv| 



M = 20 is considered. It is clearly seen at t > the per- 
turbation induced by the switch-on of the bias (Ul ^ 
and Ur — 0) propagates both inside the L lead (left- 
ward) and the normal region (rightward) with velocities 
vs ~ 2ts and ~ 2i^v respectively. At long time the 
density displays stationary oscillations due to the stabi- 
lization of the ac Josephson regime. In particular on the 
left lead the average value of Sn m -t(t) is lower with re- 
spect to the one in the right lead, since Ul > Ur. In 



Fig 13 we plot the transient behavior of the charge den- 
sity for a junction in which the (equilibrium) condition 
for perfect Andreev reflection given in Eq.( |31[ ) has been 
imposed. Remarkably we see that no appreciable density 
variation inside the lead R occurs before a dwelling time 
given by 



idwell ~ n t ar . 



(44) 
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where n = 2A/(Ur — U£) and where £ar = M/vn is 
the time needed to cross the normal chain between two 
consecutive reflections. Indeed for t < tdwcii an electron 
inside the N region undergoes n (almost) perfect An- 
dreev reflections before being transmitted through the 
right interface. The pattern of these multiple reflections 
is clearly visible in Fig |13| in which the model parameters 
are chosen in order to have n = 4 and id-well ~ 36. 



VII. SUMMARY AND CONCLUSIONS 

In this paper we have studied the dc and ac transport 
properties of a tight-binding S-N-S junction. In the dc 
case we identified three contributions to the dc Josephson 
current coming from the Andreev bound states, normal 
bound states and continuum states. The calculation of 
the latter contribution has been performed by employing 
an exact embedding procedure which consists in integrat- 
ing out the superconducting degrees of freedom and in ex- 
pressing the Nambu-Gorkov Keldysh Green's function in 
terms of the embedding self-energy. For the bound-state 
contributions we calculated the phase derivative of the 
eigenenergies of all occupied discrete states. The secular 
problem is cast in terms of an effective energy-dependent 
Hamiltonian in which the on-site and pairing potentials 
of the normal chain are renormalized via the embedding 
self-energy. The bound-state eigenenergies of chains of 
arbitrary length arc determined from a general equation 
which includes the full frequency dependence of the em- 
bedding self-energy. The limit of long-chains allows for 
further analytic manipulations and the ABS's contribu- 
tion to the total dc Josephson current is expressed in 
terms of energy-dependent phase-shifts. 

For ID superconducting leads we obtain an exact for- 
mula for the embedding self-energy at half-filling. Ex- 
plicit numerical results have been presented for short 
and long chains, and different regimes have been ana- 
lyzed. The Ishii's sawtooth behavior results from a subtle 
cancellation of highly non-linear continuum and ABS's 
contributions while the normal bound-state contribution 
vanishes. For chain hoppings smaller than half of 
the superconducting order parameter A we numerically 
observed that the dc Josephson current is entirely car- 
ried by the ABS's. This circumstance has been analyt- 
ically investigated in the limit of long chains. The con- 
dition for the occurrence of the Ishii's sawtooth behav- 
ior is expressed in terms of the microscopic parameters 
of the model. We here also point out a limitation of 
the WBL approximation, i.e. the independence of the 
current-phase relation from A. 

The ac Josephson regime was studied by applying a 
constant bias voltage across the junction and solving 
numerically the time-dependent Bogoliubov-dcGcnncs 
equations for finite leads. We used the partition-free ini- 
tial conditions for which the system is contacted and in 
equilibrium before an external driving force is switched 
on. If the leads are sufficiently long the results of the 



time propagation are the same as those of a truly infinite 
systems up to a critical time at which finite size effects 
appear.^ Such critical time is, however, large enough to 
allow for studying transient responses as well as the ac 
Josephson regime setting in after all transient effects have 
been washed out. The transient time-scale is dictated by 
the dwelling time during which an electron undergoes 
several Andreev reflections before being transmitted. By 
extracting the dc component of the ac Josephson current 
we have been able to reproduce a well-defined subgap 
structure in the current-voltage characteristics of a S-S 
junction. As expected the characteristics displays kinks 
at biases ~ 2A/n. The time-dependent approach also 
permits to perform a spectral decomposition of the ac 
current. By Fourier transforming the curve I Lit) in a 
proper time window we computed both the dissipative 
and non-dissipative components. Such procedure can be 
easily generalized to arbitrary time-dependent fields like, 
e.g., ac or pulsed biases, at the same computational cost 
and pr ovides an alternative approach to Floquet-based 
schemed 11 * 21 !. We also wish to emphasize that within the 
present approach a full microscopic description of the su- 
perconductors is provided, and hence we are able describe 
the electron dynamics not only inside the normal region, 
but also inside the leads^. This allows us to gain further 
information with respect to the WBL approximation. In 
conclusion we would like to point out the proposed time- 
dependent approach is not limited to ID electrodes and 
can be readily generalized to investigate more realistic 
superconductor-normal metal interfaces. In particular 
it would be interesting to study the case in which the 
normal region is two- dime nsional, since it has been ex- 
perimentally observe d 35 * 36 ^ and theoretically predicted^ 
that in such systems the ac Josephson current displays 
a dominant Fourier component at twice the fundamental 
Josephson frequency. 
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APPENDIX A: DERIVATION OF THE 
EMBEDDING SELF-ENERGY 



For ID leads the coupling V q is given in Eq.(12l and 
therefore the retarded embedding self-energy for lead a 
in Eq.Q reads 



= jt^^2sin 2 {q) a z g r aq iuj)a z 



(Al) 



The retarded Green's function of the uncontacted a lead 
at half-filling is given in terms of the 2x2 g-dependent 
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Bogoliubov-deGennes Hamiltonian 



where / is the Fermi distribution function. 



A a e*° 



A a e~™° 



as 



(A2) 



q = 



Introducing the eigenvectors 
/ 



(A3) 



(A4) 



APPENDIX B: CALCULATION OF THE BOUND 
STATES 



I n th is Appendix we solve the eigenvalue problem in 
Eq.(21 ). We use the following ansatJ^for the eigenstate 



amplitudes tpk(j) on the j-th site of the normal region 



of H aq with eigenvalues — ±Js 2 + A 2 and taking 



the limit A — > oo Eq.( Al ) becomes 



E»=2£ f ^ sin^) 2 Y, - 

Jo 7T f~i 



The integral can be computed analytically to yield 
m a (u> + it]) A a (uj + irj)e tx 



where 



A a (w + irj)e %Xa m a (uj + irj) 



(A5) 
(A6) 



, 4 ^A 2 - z 2 - VA 2 - z 2 + 4t| 



2t| 



A a (z) 



! 24 



with z is a complex frequency. 

The other relevant components of the Nambu self- 
energy are easily obtained starting from the retarded one: 



s<H = -/HKM-SM] 



(A8) 



Ufc(i) 

Vk(j) 



A k e ihj + B k e- ik i 



(Bl) 

where j = l,...M. Due to the symmetry of the prob- 
lem the wavefunction must be chosen so as to fulfill the 
condition 



\u k (l)\ = \u k (M)\, (B2) 
which is equivalent to |ffc(l)| = |i>fc(M)|. The above con- 
dition provides an equation for the allowed wavevectors 
k, similarly to the case of normal open chains. In the 
following we specialize to even M. The case of odd 
M is similar and does not introduce extra complica- 
tions. We first observe that due to the choice Al — An 
the effective Hamiltonian in Eq.(22) is invariant un- 
der the transformation T: Cjf — > (— 1) j cm+i- 



and 



Cji — > (— I) CM+i-j] ■ It is straightforward to realize that 
T 2 = — 1 and hence the wavefunctions obey the symme- 
try constraint v k (M) = ivu k {\) where v = ± is a parity 
index. By applying the Schrodinger equation to sites 1,2 
for spin f and to sites M, M — 1 for spin J, and exploiting 
the above symmetry constraint we obtain the following 
linear systems for the coefficients A k , B k ,C k , D k 



(e 



i k 



- ik 



t N e 



2ik 



\-A k e- l */ 2 e ikM -A k e 





t N e- 2ik -A k er ix l 2 e Lk -A k e- l xl 2 e~ lk 

q gifeM g— ikM 

-ix/2 e -ikM _ tNe ik(M-l) _ tNe -ik(M-l) 



\ 




f 1 


\ 






2t^ cos(fc) — m k 












) 




\ —iv(2tN cos(fc) — 


m k ) ) 



(B3) 



where we have chosen u k (l) = A k e lk + B k eT %k = 1. system provides the fc-dependent coefficients 
Indeed the proper normalization factor of the Andreev 
bound state wavefunction is inessential to the calcula- 
tion of the bound state energy. The solution of the above 



e m+ kM ){m k -t N e ik )-ivA k e? k , 
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B k = 



5 ik(M+l) 



e* 2 (ijv - m k e ) + ivA k e 



ikM 



D k = 



where 



1 

Jk(M+l) 
fit 



Jk-M 



ive* 2 (t N - m k e ) - A k e 



lve ^+ kM ){m k - t N e lk ) + A k e lk 



n k = t N e l ^ +kM) (e 2tk - 1) - iv~A k (e 2lk - e 2lkM ) . (B4) 



Inserting the above solution in Eq.(B2l and taking into 



account the normalization condition u k = 1 one finds the 
following equation for the allowed values of k 

= 2 [t N sin(kM) - m k smk(M- l)] 2 - t% - A 2 k 
+ t 2 N cos(2fc) + A 2 cos 2k(M - 1) 
- i/4tjvA fe cos(x/2)sin(fc)sinfc(M - 1) . (B5) 



Isolating the last term and squaring, the dependence on 
v disappears and we end up with an equation valid for 
both parities 



= t% sin 2 k(M + I) + (-l) M 2t 2 N A 2 k cos x sin 2 k 

- 2t N A 2 k sin 2 (kM) + A\sin 2 k{M- 1) 

- m 2 k [2t N sin(kM) - m k sink(M - l)] 2 , (B6) 



which coincides with Eq.(23 1. The bound states eigenen- 
ergies E^ = 2t N cosk are obtained by solving Eq.( |B6| 
numerically and retaining only the values of k such 
that \E^\ < A (Andreev bound states) and \E^\ > 
y&tg + A 2 (normal bound states). We notice that the 
wavevectors k for the ABS's are real valued while for 
normal bound states are in general complex. Indeed for 
|iiv| < \ \f^s + A 2 the energy lies below/above the 
continuum only for complex k. 
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